
"""
定义一些用于调和分析的天文参数信息以及格网参数信息
"""
from math import *
import datetime
import sys
# 添加pydggs路径
root_path = r"__file__"
if root_path not in sys.path:
    sys.path.insert(0, root_path)

from PyDGGS.PyDGGS import *

level = 6                                                                       # 格网层级
grid_type = DGGSGridType_DGGS_ISEA4H                                            # 格网类型
element_type = DGGSElementType_Cell                                             # 要素类型

TIFF_NODATA_VALUE = -9999.0                                                     # 定义TIFF文件无效值
NAN = -9999.0                                                                   # 定义NAN值,用于存储无效值
MAX_LONGTITUDE = 180                                                            # 定义经度的最大值
ONE_HOUR_TO_SECOND = 3600                                                       # 定义1小时对应的秒数
HOURS_PER_DAY = 24                                                              # 定义一天对应的秒数

wave_names = ['K2', 'N2', 'S2', 'M2', 'Q1', 'O1', 'K1', 'P1']                   # 八大分潮
H_min = 0                                                                       # 定义振幅的最小值
H_step = 0.1                                                                    # 定义振幅的等值线步长
G_min = 0                                                                       # 定义迟角的最小值
G_step = 30                                                                     # 定义迟角的等值线步长

tide_start_time = datetime.datetime(1900, 1, 1)
"""
    天文变量的初始值
    初始时间为：1900年1月1日0时
"""

s0 = 277.025
h0 = 280.19
p0 = 334.385
N0 = 259.157
ps0 = 281.221
tao_0 = h0 - s0


# 杜德森数
Doodson = {
    'M2': {
        'μ1': 2, 'μ2': 0, 'μ3': 0, 'μ4': 0, 'μ5': 0, 'μ6': 0, 'μ0': 0
    },
    'S2': {
        'μ1': 2, 'μ2': 2, 'μ3': -2, 'μ4': 0, 'μ5': 0, 'μ6': 0, 'μ0': 0
    },
    'N2': {
        'μ1': 2, 'μ2': -1, 'μ3': 0, 'μ4': 1, 'μ5': 0, 'μ6': 0, 'μ0': 0
    },
    'K2': {
        'μ1': 2, 'μ2': 2, 'μ3': 0, 'μ4': 0, 'μ5': 0, 'μ6': 0, 'μ0': 0
    },
    'K1': {
        'μ1': 1, 'μ2': 1, 'μ3': 0, 'μ4': 0, 'μ5': 0, 'μ6': 0, 'μ0': 1
    },
    'O1': {
        'μ1': 1, 'μ2': -1, 'μ3': 0, 'μ4': 0, 'μ5': 0, 'μ6': 0, 'μ0': -1
    },
    'P1': {
        'μ1': 1, 'μ2': 1, 'μ3': -2, 'μ4': 0, 'μ5': 0, 'μ6': 0, 'μ0': -1
    },
    'Q1': {
        'μ1': 1, 'μ2': -2, 'μ3': 0, 'μ4': 1, 'μ5': 0, 'μ6': 0, 'μ0': -1
    }
}

σ = {
    'M2': Doodson['M2']['μ1'] * 14.49205211 + Doodson['M2']['μ2'] * 0.54901653 + Doodson['M2']['μ3'] * 0.04106864
          + Doodson['M2']['μ4'] * 0.00464183 + Doodson['M2']['μ5'] * 0.00220641 + Doodson['M2']['μ6'] * 0.000001961,
    'N2': Doodson['N2']['μ1'] * 14.49205211 + Doodson['N2']['μ2'] * 0.54901653 + Doodson['N2']['μ3'] * 0.04106864
          + Doodson['N2']['μ4'] * 0.00464183 + Doodson['N2']['μ5'] * 0.00220641 + Doodson['N2']['μ6'] * 0.000001961,
    'S2': Doodson['S2']['μ1'] * 14.49205211 + Doodson['S2']['μ2'] * 0.54901653 + Doodson['S2']['μ3'] * 0.04106864
          + Doodson['S2']['μ4'] * 0.00464183 + Doodson['S2']['μ5'] * 0.00220641 + Doodson['S2']['μ6'] * 0.000001961,
    'K2': Doodson['K2']['μ1'] * 14.49205211 + Doodson['K2']['μ2'] * 0.54901653 + Doodson['K2']['μ3'] * 0.04106864
          + Doodson['K2']['μ4'] * 0.00464183 + Doodson['K2']['μ5'] * 0.00220641 + Doodson['K2']['μ6'] * 0.000001961,
    'K1': Doodson['K1']['μ1'] * 14.49205211 + Doodson['K1']['μ2'] * 0.54901653 + Doodson['K1']['μ3'] * 0.04106864
          + Doodson['K1']['μ4'] * 0.00464183 + Doodson['K1']['μ5'] * 0.00220641 + Doodson['K1']['μ6'] * 0.000001961,
    'O1': Doodson['O1']['μ1'] * 14.49205211 + Doodson['O1']['μ2'] * 0.54901653 + Doodson['O1']['μ3'] * 0.04106864
          + Doodson['O1']['μ4'] * 0.00464183 + Doodson['O1']['μ5'] * 0.00220641 + Doodson['O1']['μ6'] * 0.000001961,
    'P1': Doodson['P1']['μ1'] * 14.49205211 + Doodson['P1']['μ2'] * 0.54901653 + Doodson['P1']['μ3'] * 0.04106864
          + Doodson['P1']['μ4'] * 0.00464183 + Doodson['P1']['μ5'] * 0.00220641 + Doodson['P1']['μ6'] * 0.000001961,
    'Q1': Doodson['Q1']['μ1'] * 14.49205211 + Doodson['Q1']['μ2'] * 0.54901653 + Doodson['Q1']['μ3'] * 0.04106864
          + Doodson['Q1']['μ4'] * 0.00464183 + Doodson['Q1']['μ5'] * 0.00220641 + Doodson['Q1']['μ6'] * 0.000001961,
}
'''
    σ 为各分潮的角速度
    计算公式为：σ = μ1*σ_tao + μ2*σ_s + μ3*σ_h + μ4*σ_p + μ5*σ_n + μ6*σ_ps
'''

t1 = datetime.datetime.strptime('1950-1-1 00:00:00', '%Y-%m-%d %H:%M:%S')
t2 = datetime.datetime.strptime('1900-1-1 00:00:00', '%Y-%m-%d %H:%M:%S')

delta = (t1 - t2).days

V0 = {
    'M2': Doodson['M2']['μ1'] * tao_0 + Doodson['M2']['μ2'] * s0 + Doodson['M2']['μ3'] * h0 + Doodson['M2']['μ4'] * p0 +
          Doodson['M2']['μ5'] * N0 + Doodson['M2']['μ6'] * ps0 + Doodson['M2']['μ0'] * 90,
    'S2': Doodson['S2']['μ1'] * tao_0 + Doodson['S2']['μ2'] * s0 + Doodson['S2']['μ3'] * h0 + Doodson['S2']['μ4'] * p0 +
          Doodson['S2']['μ5'] * N0 + Doodson['S2']['μ6'] * ps0 + Doodson['S2']['μ0'] * 90,
    'N2': Doodson['N2']['μ1'] * tao_0 + Doodson['N2']['μ2'] * s0 + Doodson['N2']['μ3'] * h0 + Doodson['N2']['μ4'] * p0 +
          Doodson['N2']['μ5'] * N0 + Doodson['N2']['μ6'] * ps0 + Doodson['N2']['μ0'] * 90,
    'K2': Doodson['K2']['μ1'] * tao_0 + Doodson['K2']['μ2'] * s0 + Doodson['K2']['μ3'] * h0 + Doodson['K2']['μ4'] * p0 +
          Doodson['K2']['μ5'] * N0 + Doodson['K2']['μ6'] * ps0 + Doodson['K2']['μ0'] * 90,
    'K1': Doodson['K1']['μ1'] * tao_0 + Doodson['K1']['μ2'] * s0 + Doodson['K1']['μ3'] * h0 + Doodson['K1']['μ4'] * p0 +
          Doodson['K1']['μ5'] * N0 + Doodson['K1']['μ6'] * ps0 + Doodson['K1']['μ0'] * 90,
    'O1': Doodson['O1']['μ1'] * tao_0 + Doodson['O1']['μ2'] * s0 + Doodson['O1']['μ3'] * h0 + Doodson['O1']['μ4'] * p0 +
          Doodson['O1']['μ5'] * N0 + Doodson['O1']['μ6'] * ps0 + + Doodson['Q1']['μ0'] * 90,
    'P1': Doodson['P1']['μ1'] * tao_0 + Doodson['P1']['μ2'] * s0 + Doodson['P1']['μ3'] * h0 + Doodson['P1']['μ4'] * p0 +
          Doodson['P1']['μ5'] * N0 + Doodson['P1']['μ6'] * ps0 + + Doodson['P1']['μ0'] * 90,
    'Q1': Doodson['Q1']['μ1'] * tao_0 + Doodson['Q1']['μ2'] * s0 + Doodson['Q1']['μ3'] * h0 + Doodson['Q1']['μ4'] * p0 +
          Doodson['Q1']['μ5'] * N0 + Doodson['Q1']['μ6'] * ps0 + + Doodson['Q1']['μ0'] * 90,
    }
'''
    V0:天文初相角
    计算公式：V0 = μ1*tao_0 + μ2*s_0 + μ3*h_0 + μ4*p_0 + μ5*n_0 + μ6*ps_0 + μ0*90
    初始时间为:1900年1月1日0时
'''

def cal_v(t):
    """
        任意时刻的天文相角V[分潮] = σ[分潮] * t + V0
    """
    v = {
        'M2': σ['M2'] * t + V0['M2'],
        'S2': σ['S2'] * t + V0['S2'],
        'N2': σ['N2'] * t + V0['N2'],
        'K2': σ['K2'] * t + V0['K2'],
        'K1': σ['K1'] * t + V0['K1'],
        'O1': σ['O1'] * t + V0['O1'],
        'P1': σ['P1'] * t + V0['P1'],
        'Q1': σ['Q1'] * t + V0['Q1'],
    }

    return v


def cal_f_u(t):
    """
        计算订正因子f、订正角u
        计算方法：采用陈宗镛发表于1900年《潮汐分析和推算的一种j、u模型》
    """
    p = radians(0.00464183 * t + p0)
    n = radians(-0.00220641 * t + N0)
    ps = radians(0.00000196 * t + ps0)

    fu_czy = {
    "M2": {
        'fcosu': 1 - 0.03733 * cos(n) + 0.0052 * cos(2 * n) + 0.0058 * cos(2 * p) + 0.00021 * cos(2 * p - n),
        'fsinu': -0.03733 * sin(n) + 0.0052 * sin(2 * n) + 0.0058 * sin(2 * p) + 0.00021 * sin(2 * p - n)
    },
    "S2": {
        'fcosu': 1 + 0.00225 * cos(n) + 0.00014 * cos(2 * p),
        'fsinu': 0.00225 * sin(n) + 0.00014 * sin(2 * p)
    },
    "N2": {
        'fcosu': 1 - 0.03733 * cos(n) + 0.00052 * cos(2 * n) + 0.00081 * cos(p - ps) - 0.00385 * cos(2 * p - 2 * n),
        'fsinu': -0.03733 * sin(n) + 0.00052 * sin(2 * n) - 0.00081 * sin(p - ps) + 0.00365 * sin(2 * p - 2 * n)
    },
    "K2": {
        'fcosu': 1 + 0.28518 * cos(n) + 0.03235 * cos(2 * n),
        'fsinu': -0.31074 * sin(n) - 0.03235 * sin(2 * n)
    },
    "K1": {
        'fcosu': 1 + 0.11573 * cos(n) - 0.00281 * cos(2 * n) + 0.00019 * cos(2 * p - n),
        'fsinu': -0.15539 * sin(n) + 0.00303 * sin(2 * n) - 0.00019 * sin(2 * p - n)
    },
    "O1": {
        'fcosu': 1 + 0.18852 * cos(n) - 0.00578 * cos(2 * n) - 0.00645 * cos(2 * p) - 0.00103 * cos(2 * p - n) + 0.00019 * cos(2 * p + n),
        'fsinu': 0.18852 * sin(n) - 0.00578 * sin(2 * n) - 0.00645 * sin(2 * p) - 0.00103 * sin(2 * p - n) + 0.00019 * sin(2 * p + n)
    },
    "P1": {
        'fcosu': 1 - 0.01123 * cos(n) - 0.00040 * cos(2 * ps) - 0.00148 * cos(2 * p) - 0.00029 * cos(2 * p - n) + 0.00080 * cos(2 * n),
        'fsinu': - 0.01123 * sin(n) - 0.00040 * sin(2 * ps) - 0.00148 * sin(2 * p) - 0.00029 * sin(2 * p - n) + 0.00080 * sin(2 * n)
    },
    "Q1": {
        'fcosu': 1 + 0.18844 * cos(n) - 0.00568 * cos(2 * n) - 0.00277 * cos(2 * p) - 0.00388 * cos(2 * p - 2*n) + 0.00083 * cos(p - ps) - 0.00069 *cos(2*p -3*n),
        'fsinu': 0.18844 * sin(n) - 0.00568 * sin(2 * n) - 0.00277 * sin(2 * p) + 0.00388 * sin(2 * p - 2*n) - 0.00083 * sin(p - ps) + 0.00069 *sin(2*p -3*n)
    }

    }

    return fu_czy
